inter.binning.90 <- function (Y, D, X, Z = NULL, FE = NULL, weights = NULL, data, 
          na.rm = FALSE, nbins = 3, cutoffs = NULL, vartype = "homoscedastic", 
          cl = NULL, time = NULL, pairwise = TRUE, figure = TRUE, main = NULL, 
          Ylabel = NULL, Dlabel = NULL, Xlabel = NULL, xlim = NULL, 
          ylim = NULL, interval = NULL, Xdistr = "histogram", wald = TRUE) 
{
  x <- NULL
  y <- NULL
  xmin <- NULL
  xmax <- NULL
  ymin <- NULL
  ymax <- NULL
  if (is.character(Y) == FALSE) {
    stop("Y is not a string.")
  }
  else {
    Y <- Y[1]
  }
  if (is.character(D) == FALSE) {
    stop("D is not a string.")
  }
  else {
    D <- D[1]
  }
  if (is.character(X) == FALSE) {
    stop("X is not a string.")
  }
  else {
    X <- X[1]
  }
  if (is.null(Z) == FALSE) {
    for (i in 1:length(Z)) {
      if (is.character(Z[i]) == FALSE) {
        stop("Some element in Z is not a string.")
      }
    }
  }
  if (is.null(FE) == FALSE) {
    for (i in 1:length(FE)) {
      if (is.character(FE[i]) == FALSE) {
        stop("Some element in FE is not a string.")
      }
    }
  }
  if (is.null(weights) == FALSE) {
    if (is.character(weights) == FALSE) {
      stop("weigths is not a string.")
    }
    else {
      weights <- weights[1]
    }
  }
  if (is.data.frame(data) == FALSE) {
    stop("Not a data frame.")
  }
  if (is.logical(na.rm) == FALSE & is.numeric(na.rm) == FALSE) {
    stop("na.rm is not a logical flag.")
  }
  if (is.null(nbins) == FALSE) {
    if (nbins%%1 != 0) {
      stop("nbins is not a positive integer.")
    }
    else {
      nbins <- nbins[1]
    }
    if (nbins < 1) {
      stop("nbins is not a positive integer.")
    }
  }
  if (is.null(cutoffs) == FALSE) {
    if (is.numeric(cutoffs) == FALSE) {
      stop("Some element of cutoffs is not numeric.")
    }
  }
  if (is.null(vartype) == TRUE) {
    vartype <- "homoscedastic"
  }
  if (!vartype %in% c("homoscedastic", "robust", "cluster", 
                      "pcse")) {
    stop("vartype must be one of the following")
  }
  else if (vartype == "cluster") {
    if (is.null(cl) == TRUE) {
      stop("cl not specified")
    }
  }
  else if (vartype == "pcse") {
    if (is.null(cl) == TRUE | is.null(time) == TRUE) {
      stop("cl or time not specified; set cl =")
    }
  }
  if (is.null(cl) == FALSE) {
    if (is.character(cl) == FALSE) {
      stop("cl is not a string.")
    }
    else {
      cl <- cl[1]
      if (vartype != "pcse") {
        vartype <- "cluster"
      }
    }
  }
  if (is.null(time) == FALSE) {
    if (is.character(time) == FALSE) {
      stop("time is not a string.")
    }
    else {
      time <- time[1]
    }
  }
  if (is.logical(pairwise) == FALSE & is.numeric(pairwise) == 
      FALSE) {
    stop("pairwise is not a logical flag.")
  }
  if (is.null(Ylabel) == TRUE) {
    Ylabel <- Y
  }
  else {
    if (is.character(Ylabel) == FALSE) {
      stop("Ylabel is not a string.")
    }
    else {
      Ylabel <- Ylabel[1]
    }
  }
  if (is.null(Dlabel) == TRUE) {
    Dlabel <- D
  }
  else {
    if (is.character(Dlabel) == FALSE) {
      stop("Dlabel is not a string.")
    }
    else {
      Dlabel <- Dlabel[1]
    }
  }
  if (is.null(Xlabel) == TRUE) {
    Xlabel <- X
  }
  else {
    if (is.character(Xlabel) == FALSE) {
      stop("Xlabel is not a string.")
    }
    else {
      Xlabel <- Xlabel[1]
    }
  }
  if (is.null(main) == FALSE) {
    main <- main[1]
  }
  if (!Xdistr %in% c("hist", "histogram", "density")) {
    stop("Xdistr must be either")
  }
  if (is.null(xlim) == FALSE) {
    if (is.numeric(xlim) == FALSE) {
      stop("Some element in xlim is not numeric.")
    }
    else {
      if (length(xlim) != 2) {
        stop("xlim must be of length 2.")
      }
    }
  }
  if (is.null(ylim) == FALSE) {
    if (is.numeric(ylim) == FALSE) {
      stop("Some element in ylim is not numeric.")
    }
    else {
      if (length(ylim) != 2) {
        stop("ylim must be of length 2.")
      }
    }
  }
  if (na.rm == TRUE) {
    data <- na.omit(data[, c(Y, D, X, Z, FE)])
  }
  else {
    if (sum(is.na(data[, c(Y, D, X, Z, FE)])) > 0) {
      stop("Missing values. Try option na.rm = TRUE")
    }
  }
  requireNamespace("ggplot2")
  n <- dim(data)[1]
  data[, D] <- as.numeric(data[, D])
  if (is.null(FE) == FALSE) {
    if (is.null(Z) == TRUE) {
      Z <- c()
    }
    for (i in 1:length(FE)) {
      Z <- c(Z, paste("as.factor(", FE[i], ")", sep = ""))
    }
  }
  if (is.null(Z) == FALSE) {
    mod.f <- as.formula(paste(Y, "~", D, "+", X, "+", D, 
                              "*", X, "+", paste(Z, collapse = "+"), sep = ""))
  }
  else {
    mod.f <- as.formula(paste(Y, "~", D, "+", X, "+", D, 
                              "*", X, sep = ""))
  }
  if (is.null(weights) == TRUE) {
    mod.naive <- lm(mod.f, data = data)
  }
  else {
    mod.naive <- lm(mod.f, data = data, weights = data[, 
                                                       weights])
  }
  coefs <- summary(mod.naive)$coefficients[, 1]
  coef.D <- coefs[D]
  coef.X <- coefs[X]
  coef.DX <- coefs[paste(D, X, sep = ":")]
  if (is.null(vartype) == TRUE) {
    vartype <- "homoscedastic"
  }
  if (vartype == "homoscedastic") {
    v <- vcov(mod.naive)
  }
  else if (vartype == "robust") {
    requireNamespace("sandwich")
    v <- vcov(mod.naive, type = "HC1")
  }
  else if (vartype == "cluster") {
    v <- vcovCluster(mod.naive, cluster = data[, cl])
  }
  else if (vartype == "pcse") {
    requireNamespace("pcse")
    v <- pcse(mod.naive, groupN = data[, cl], groupT = data[, 
                                                            time], pairwise = pairwise)$vcov
  }
  if (vartype == "pcse") {
    var.D <- v[D, D]
    var.DX <- v[paste(D, X, sep = "."), paste(D, X, sep = ".")]
    cov.DX <- v[D, paste(D, X, sep = ".")]
  }
  else {
    var.D <- v[D, D]
    var.DX <- v[paste(D, X, sep = ":"), paste(D, X, sep = ":")]
    cov.DX <- v[D, paste(D, X, sep = ":")]
  }
  X.lvls <- as.numeric(quantile(data[, X], probs = seq(0, 1, 
                                                       0.01)))
  marg <- coef.D + coef.DX * X.lvls
  marg
  se <- sqrt(var.D + X.lvls^2 * var.DX + 2 * X.lvls * cov.DX)
  df <- mod.naive$df.residual
  crit <- abs(qt(0.05, df = df))
  lb <- marg - crit * se
  ub <- marg + crit * se
  if (is.null(cutoffs) == TRUE) {
    cuts.X <- quantile(data[, X], probs = seq(0, 1, 1/nbins))
    while (length(unique(cuts.X)) != nbins + 1) {
      nbins <- nbins - 1
      cuts.X <- quantile(data[, X], probs = seq(0, 1, 1/nbins))
    }
  }
  else {
    cutoffs <- cutoffs[which(cutoffs > min(data[, X]) & cutoffs < 
                               max(data[, X]))]
    cuts.X <- sort(unique(c(min(data[, X]), cutoffs, max(data[, 
                                                              X]))))
  }
  groupX <- cut(data[, X], breaks = cuts.X, labels = FALSE)
  groupX[which(data[, X] == min(data[, X]))] <- 1
  nbins <- length(unique(groupX))
  groupX2 <- cut(data[, X], breaks = cuts.X)
  gp.lab = paste(Xlabel, ": ", levels(groupX2), sep = "")
  gp.lab[1] <- paste(Xlabel, ": [", substring(levels(groupX2)[1], 
                                              2), sep = "")
  x0 <- rep(NA, nbins)
  for (i in 1:nbins) x0[i] <- median(data[which(groupX == i), 
                                          X], na.rm = TRUE)
  G <- DG <- GX <- DGX <- matrix(0, n, nbins)
  for (i in 1:nbins) {
    G[which(groupX == i), i] <- 1
    DG[, i] <- data[, D] * G[, i]
    GX[, i] <- G[, i] * (data[, X] - x0[i])
    DGX[, i] <- DG[, i] * (data[, X] - x0[i])
  }
  Gs <- GXs <- DGs <- DGXs <- c()
  for (i in 1:nbins) {
    Gs <- c(Gs, paste("G[,", i, "]", sep = ""))
    GXs <- c(GXs, paste("GX[,", i, "]", sep = ""))
    DGs <- c(DGs, paste("DG[,", i, "]", sep = ""))
    DGXs <- c(DGXs, paste("DGX[,", i, "]", sep = ""))
  }
  Xf <- paste(Y, "~ -1+", paste(DGs, collapse = "+"), "+", 
              paste(DGXs, collapse = "+"), "+", paste(Gs, collapse = "+"), 
              "+", paste(GXs, collapse = "+"), sep = "")
  if (is.null(Z) == FALSE) {
    Xf <- paste(Xf, "+", paste(Z, collapse = "+"), sep = "")
  }
  mod.Xf <- as.formula(Xf)
  if (is.null(weights) == TRUE) {
    mod.X <- lm(mod.Xf, data = data)
  }
  else {
    mod.X <- lm(mod.Xf, data = data, weights = data[, weights])
  }
  Xcoefs <- mod.X$coefficients[1:nbins]
  if (vartype == "homoscedastic") {
    X.v <- vcov(mod.X)
  }
  else if (vartype == "robust") {
    X.v <- vcov(mod.X, type = "HC1")
  }
  else if (vartype == "cluster") {
    X.v <- vcovCluster(mod.X, cluster = data[, cl])
  }
  else if (vartype == "pcse") {
    if (is.null(Z) == FALSE) {
      exclude <- names(which(is.na(mod.X$coefficients) == 
                               TRUE))
      Z.ex <- setdiff(Z, exclude)
      Xf <- paste(Y, "~ -1+", paste(DGs, collapse = "+"), 
                  "+", paste(DGXs, collapse = "+"), "+", paste(Gs, 
                                                               collapse = "+"), "+", paste(GXs, collapse = "+"), 
                  "+", paste(Z.ex, collapse = "+"), sep = "")
      mod.X <- lm(as.formula(Xf), data = data)
    }
    X.v <- pcse(mod.X, groupN = data[, cl], groupT = data[, 
                                                          time], pairwise = pairwise)$vcov
  }
  X.v <- X.v[1:nbins, 1:nbins]
  X.se <- sqrt(diag(as.matrix(X.v, drop = FALSE)))
  X.se[which(is.na(Xcoefs))] <- NA
  df.X <- mod.X$df.residual
  crit.X <- abs(qt(0.05, df = df.X))
  lb.X <- Xcoefs - crit.X * X.se
  ub.X <- Xcoefs + crit.X * X.se
  est.binning <- cbind(x0, coef = Xcoefs, se = X.se, CI_lower = lb.X, 
                       CI_upper = ub.X)
  rownames(est.binning) <- gp.lab
  if (figure == TRUE) {
    if (is.null(Xlabel) == FALSE) {
      x.label <- c(paste("Moderator: ", Xlabel, sep = ""))
      y.label <- c(paste("Marginal Effect of Treatment"))
    }
    else {
      x.label <- c(paste("Moderator: ", X, sep = ""))
      y.label <- c(paste("Marginal Effect of Treatment"))
    }
    out <- data.frame(X.lvls, marg, lb, ub)
    out.bin <- data.frame(x0, Xcoefs, lb.X, ub.X)
    out.bin2 <- out.bin[which(is.na(Xcoefs) == FALSE), ]
    out.bin3 <- out.bin[which(is.na(Xcoefs) == TRUE), ]
    errorbar.width <- (max(X.lvls) - min(X.lvls))/20
    p1 <- ggplot()
    yrange <- na.omit(c(marg, lb, ub, Xcoefs, lb.X, ub.X))
    if (is.null(ylim) == FALSE) {
      yrange <- c(ylim[2], ylim[1] + (ylim[2] - ylim[1]) * 
                    1/6)
    }
    maxdiff <- (max(yrange) - min(yrange))
    pos <- max(yrange) - maxdiff/20
    p1 <- p1 + geom_hline(yintercept = 0, colour = "white", 
                          size = 2)
    if (is.null(interval) == FALSE) {
      p1 <- p1 + geom_vline(xintercept = interval, colour = "steelblue", 
                            linetype = 2, size = 1.5)
    }
    if (is.null(Xdistr) == TRUE) {
      Xdistr <- "density"
    }
    else if (Xdistr != "density" & Xdistr != "histogram" & 
             Xdistr != "hist") {
      Xdistr <- "density"
    }
    if (Xdistr == "density") {
      if (length(unique(data[, D])) == 2) {
        if (is.null(weights) == TRUE) {
          de.co <- density(data[data[, D] == 0, X])
          de.tr <- density(data[data[, D] == 1, X])
        }
        else {
          suppressWarnings(de.co <- density(data[data[, 
                                                      D] == 0, X], weights = data[data[, D] == 
                                                                                    0, weights]))
          suppressWarnings(de.tr <- density(data[data[, 
                                                      D] == 1, X], weights = data[data[, D] == 
                                                                                    1, weights]))
        }
        deX.ymin <- min(yrange) - maxdiff/5
        deX.co <- data.frame(x = de.co$x, y = de.co$y/max(de.co$y) * 
                               maxdiff/5 + min(yrange) - maxdiff/5)
        deX.tr <- data.frame(x = de.tr$x, y = de.tr$y/max(de.tr$y) * 
                               maxdiff/5 + min(yrange) - maxdiff/5)
        feed.col <- col2rgb("gray50")
        col.co <- rgb(feed.col[1]/1000, feed.col[2]/1000, 
                      feed.col[3]/1000)
        col.tr <- rgb(red = 1, blue = 0, green = 0)
        p1 <- p1 + geom_ribbon(data = deX.co, aes(x = x, 
                                                  ymax = y, ymin = deX.ymin), fill = col.co, 
                               alpha = 0.2) + geom_ribbon(data = deX.tr, aes(x = x, 
                                                                             ymax = y, ymin = deX.ymin), fill = col.tr, 
                                                          alpha = 0.2)
      }
      else {
        if (is.null(weights) == TRUE) {
          de <- density(data[, X])
        }
        else {
          suppressWarnings(de <- density(data[, X], weights = data[, 
                                                                   weights]))
        }
        deX.ymin <- min(yrange) - maxdiff/5
        deX <- data.frame(x = de$x, y = de$y/max(de$y) * 
                            maxdiff/5 + min(yrange) - maxdiff/5)
        feed.col <- col2rgb("gray50")
        col <- rgb(feed.col[1]/1000, feed.col[2]/1000, 
                   feed.col[3]/1000)
        p1 <- p1 + geom_ribbon(data = deX, aes(x = x, 
                                               ymax = y, ymin = deX.ymin), fill = col, alpha = 0.2)
      }
    }
    else {
      if (length(unique(data[, D])) == 2) {
        if (is.null(weights) == TRUE) {
          hist.out <- hist(data[, X], breaks = 80, plot = FALSE)
        }
        else {
          suppressWarnings(hist.out <- hist(data[, X], 
                                            data[, weights], breaks = 80, plot = FALSE))
        }
        n.hist <- length(hist.out$mids)
        dist <- hist.out$mids[2] - hist.out$mids[1]
        hist.max <- max(hist.out$counts)
        count1 <- rep(0, n.hist)
        treat <- which(data[, D] == max(data[, D]))
        for (i in 1:n.hist) {
          count1[i] <- sum(data[treat, X] >= hist.out$breaks[i] & 
                             data[treat, X] < hist.out$breaks[(i + 1)])
        }
        count1[n.hist] <- sum(data[treat, X] >= hist.out$breaks[n.hist] & 
                                data[treat, X] <= hist.out$breaks[n.hist + 
                                                                    1])
        histX <- data.frame(ymin = rep(min(yrange) - 
                                         maxdiff/5, n.hist), ymax = hist.out$counts/hist.max * 
                              maxdiff/5 + min(yrange) - maxdiff/5, xmin = hist.out$mids - 
                              dist/2, xmax = hist.out$mids + dist/2, count1 = count1/hist.max * 
                              maxdiff/5 + min(yrange) - maxdiff/5)
        p1 <- p1 + geom_rect(data = histX, aes(xmin = xmin, 
                                               xmax = xmax, ymin = ymin, ymax = ymax), colour = "gray50", 
                             alpha = 0, size = 0.5) + geom_rect(data = histX, 
                                                                aes(xmin = xmin, xmax = xmax, ymin = ymin, 
                                                                    ymax = count1), fill = "red", colour = "grey50", 
                                                                alpha = 0.3, size = 0.5)
      }
      else {
        hist.out <- hist(data[, X], breaks = 80, plot = FALSE)
        n.hist <- length(hist.out$mids)
        dist <- hist.out$mids[2] - hist.out$mids[1]
        hist.max <- max(hist.out$counts)
        histX <- data.frame(ymin = rep(min(yrange) - 
                                         maxdiff/5, n.hist), ymax = hist.out$counts/hist.max * 
                              maxdiff/5 + min(yrange) - maxdiff/5, xmin = hist.out$mids - 
                              dist/2, xmax = hist.out$mids + dist/2)
        p1 <- p1 + geom_rect(data = histX, aes(xmin = xmin, 
                                               xmax = xmax, ymin = ymin, ymax = ymax), colour = "gray50", 
                             alpha = 0, size = 0.5)
      }
    }
    if (is.null(main) == FALSE) {
      p1 <- p1 + ggtitle(main) + theme(plot.title = element_text(hjust = 0.5, 
                                                                 size = 12, lineheight = 0.8, face = "bold"))
    }
    if (nbins == 3) {
      p1 <- p1 + annotate(geom = "text", x = out.bin[1, 
                                                     1], y = pos, label = "L", colour = "gray50", 
                          size = 10) + annotate(geom = "text", x = out.bin[2, 
                                                                           1], y = pos, label = "M", colour = "gray50", 
                                                size = 10) + annotate(geom = "text", x = out.bin[3, 
                                                                                                 1], y = pos, label = "H", colour = "gray50", 
                                                                      size = 10)
    }
    else if (nbins == 4) {
      p1 <- p1 + annotate(geom = "text", x = out.bin[1, 
                                                     1], y = pos, label = "Q1", colour = "gray50", 
                          size = 4) + annotate(geom = "text", x = out.bin[2, 
                                                                           1], y = pos, label = "Q2", colour = "gray50", 
                                                size = 4) + annotate(geom = "text", x = out.bin[3, 
                                                                                                 1], y = pos, label = "Q3", colour = "gray50", 
                                                                      size = 4) + annotate(geom = "text", x = out.bin[4, 
                                                                                                                       1], y = pos, label = "Q4", colour = "gray50", 
                                                                                            size = 4)
    }
    p1 <- p1 + geom_line(data = out, aes(X.lvls, marg)) + 
      geom_ribbon(data = out, aes(x = X.lvls, ymin = lb, 
                                  ymax = ub), alpha = 0.2) + xlab(x.label) + ylab(y.label)
    p1 <- p1 + geom_errorbar(data = out.bin2, aes(x = x0, 
                                                  ymin = lb.X, ymax = ub.X), colour = "red", width = errorbar.width) + 
      geom_point(data = out.bin2, aes(x0, Xcoefs), size = 4, 
                 shape = 21, fill = "white", colour = "red") + 
      theme(axis.title = element_text(size = 11))
    p1 <- p1 + annotate(geom = "text", x = out.bin3[, 1], 
                        y = rep(0, dim(out.bin3)[1]), label = "NaN", colour = "red")
    if (is.null(ylim) == FALSE) {
      ylim2 = c(ylim[1] - (ylim[2] - ylim[1]) * 0.25/6, 
                ylim[2] + (ylim[2] - ylim[1]) * 0.4/6)
    }
    if (is.null(xlim) == FALSE & is.null(ylim) == FALSE) {
      p1 <- p1 + coord_cartesian(xlim = xlim, ylim = ylim2)
    }
    if (is.null(xlim) == TRUE & is.null(ylim) == FALSE) {
      p1 <- p1 + coord_cartesian(ylim = ylim2)
    }
    if (is.null(xlim) == FALSE & is.null(ylim) == TRUE) {
      p1 <- p1 + coord_cartesian(xlim = xlim)
    }
  }
  btreat <- length(unique(data[, D])) == 2
  varD <- c()
  for (i in 1:nbins) {
    varD <- c(varD, var(data[groupX == i, D]/mean(data[groupX == 
                                                         i, D])))
  }
  requireNamespace("Lmoments")
  Lkurtosis <- Lmoments(data[, X], returnobject = TRUE)$ratios[4]
  correctOrder <- ifelse(as.numeric((Xcoefs[1] - Xcoefs[2]) * 
                                      (Xcoefs[2] - Xcoefs[3])) > 0, TRUE, FALSE)
  pvalue <- function(i, j) {
    stat <- (Xcoefs[i] - Xcoefs[j])/sqrt(X.v[i, i] + X.v[j, 
                                                         j] - 2 * X.v[i, j])
    p <- (1 - pt(abs(stat), df.X)) * 2
    return(p)
  }
  if (nbins == 3) {
    p.twosided <- round(c(pvalue(1, 2), pvalue(2, 3), pvalue(1, 
                                                             3)), digits = 4)
    names(p.twosided) <- c("p.1v2", "p.2v3", "p.1v3")
    names(Xcoefs) <- c("X_low", "X_med", "X_high")
  }
  else if (nbins == 2) {
    p.twosided <- round(pvalue(1, 2), digits = 4)
    names(p.twosided) <- c("p.LvH")
    names(Xcoefs) <- c("X_low", "X_high")
  }
  else if (nbins == 4) {
    names(Xcoefs) <- c("X_low", "X_med1", "X_med2", "X_high")
  }
  formula0 <- paste(Y, "~", D, "+", X, "+", D, "*", X)
  G <- DG <- GX <- DGX <- matrix(0, n, (nbins - 1))
  for (i in 1:(nbins - 1)) {
    G[which(groupX == (i + 1)), i] <- 1
    DG[, i] <- data[, D] * G[, i]
    GX[, i] <- data[, X] * G[, i]
    DGX[, i] <- data[, D] * data[, X] * G[, i]
  }
  Gs <- GXs <- DGs <- DGXs <- c()
  for (i in 2:nbins) {
    Gs <- c(Gs, paste("G", i, sep = ""))
    GXs <- c(GXs, paste("GX", i, sep = ""))
    DGs <- c(DGs, paste("DG", i, sep = ""))
    DGXs <- c(DGXs, paste("DGX", i, sep = ""))
  }
  colnames(G) <- Gs
  colnames(DG) <- DGs
  colnames(GX) <- GXs
  colnames(DGX) <- DGXs
  data.aug <- cbind.data.frame(data, G, DG, GX, DGX)
  formula1 <- paste(formula0, "+", paste(Gs, collapse = " + "), 
                    "+", paste(GXs, collapse = " + "), "+", paste(DGs, collapse = " + "), 
                    "+", paste(DGXs, collapse = " + "))
  if (is.null(Z) == FALSE) {
    formula0 <- paste(formula0, "+", paste(Z, collapse = " + "))
    formula1 <- paste(formula1, "+", paste(Z, collapse = " + "))
  }
  if (is.null(weights) == TRUE) {
    mod.re <- lm(as.formula(formula0), data = data.aug)
    mod.un <- lm(as.formula(formula1), data = data.aug)
  }
  else {
    mod.re <- lm(as.formula(formula0), data = data.aug, weights = data.aug[, 
                                                                           weights])
    mod.un <- lm(as.formula(formula1), data = data.aug, weights = data.aug[, 
                                                                           weights])
  }
  if (is.null(vartype) == TRUE) 
    vartype <- "homoscedastic"
  if (vartype == "homoscedastic") {
    v <- vcov(mod.un)
  }
  else if (vartype == "robust") {
    v <- vcov(mod.un, type = "HC1")
  }
  else if (vartype == "cluster") {
    v <- vcovCluster(mod.un, cluster = data.aug[, cl])
  }
  else if (vartype == "pcse") {
    v <- pcse(mod.un, groupN = data.aug[, cl], groupT = data.aug[, 
                                                                 time], pairwise = pairwise)$vcov
  }
  if (wald == TRUE) {
    requireNamespace("lmtest")
    wald.out <- tryCatch(p.wald <- round(waldtest(mod.re, 
                                                  mod.un, test = "Chisq", vcov = v)[[4]][2], 4), error = function(e) {
                                                    return(NULL)
                                                  })
    if (is.null(wald.out) == TRUE) {
      p.wald <- NULL
      warning("Var-cov matrix nearly sigular in the Wald test.")
    }
  }
  out <- list(est.binning = round(est.binning, 4), binary.treatment = round(btreat, 
                                                                            4), bin.size = round(c(table(groupX)/length(groupX)), 
                                                                                                 4), treat.variation.byGroup = round(varD, 4), X.Lkurtosis = round(Lkurtosis, 
                                                                                                                                                                   4))
  if (nbins == 3) {
    out <- c(out, list(correctOrder = correctOrder))
  }
  if (nbins %in% c(2, 3)) {
    out <- c(out, list(p.twosided = p.twosided))
  }
  if (wald == TRUE) {
    out <- c(out, list(p.wald = p.wald))
  }
  if (figure == TRUE) {
    out <- c(out, list(graph = p1))
  }
  return(out)
}